Matching

Data

We’ll load up the Lalonde data. Our goal here is to see if we can replicate or at least get close to the correct estimate in Lalonde’s 1986 paper using observational data. According to his model, the actual effect for the job training program was a $1,794 increase in the subject’s wages in 1978.

As usual a good starting point here is to get some basic descriptive stats on our variables.

lalonde|>
  select(treat, educ, married, nodegree, race)|>
  mutate(treat = factor(treat, labels=c("control", "treatment")),
         married = factor(married, labels=c('single','married')),
         nodegree = factor(nodegree, labels=c('degree','no degree'))
         )|>
  ggpairs(aes(fill=treat, color=treat), upper='blank') +
  theme_bw() +
  scale_fill_brewer(palette='Dark2') +
  scale_color_brewer(palette='Dark2')

lalonde|>
  select(treat, age, re74, re75, re78)|>
  mutate(treat = factor(treat, labels=c("control", "treatment")))|>
  ggpairs(aes(fill=treat, color=treat), upper='blank') +
  theme_bw() +
  scale_fill_brewer(palette='Dark2') +
  scale_color_brewer(palette='Dark2')

Regression

We’ll start with regression based estimates. The first model includes no controls, and the second includes controls for age, education, marital status, past income, degree status, and race:

# the uncontrolled regression
mod0<-lm(re78 ~ treat, data=lalonde)
mod1<-lm(re78 ~ treat + age + educ + married  + re74 + re75+ nodegree +race,data=lalonde)


modelsummary(list("treatment only"  =mod0, 
                        "treatment + controls" = mod1),
             estimate  = "{estimate}",  
             fmt = fmt_significant(),
             statistic ='conf.int',
             conf_level = .95,
             gof_omit = 'F|RMSE|R2$|AIC|Log.Lik.', # remove some model stats 
             # add a title
             title = "DV: Wages"
             )
tinytable_4qpkbe7ylke71xo38ci0
DV: Wages
treatment only treatment + controls
(Intercept) 6984 -1174
[6276, 7693] [-5998, 3650]
treat -635 1548.2
[-1926, 655] [13.9, 3082.6]
age 13.0
[-50.8, 76.8]
educ 403.9
[91.9, 716.0]
married 407
[-959, 1772]
re74 0.2964
[0.1819, 0.4108]
re75 0.2315
[0.0261, 0.4370]
nodegree 260
[-1404, 1924]
racehispan 1740
[-261, 3740]
racewhite 1241
[-269, 2750]
Num.Obs. 614 614
R2 Adj. 0.000 0.135
BIC 12712.0 12666.1

Matching methods

Mahalanobis Scores

You can get Mahalanobis score matches by using method="nearest" and distance='mahalanobis'

ma_matched <- matchit(treat ~ age + educ +  married + race +
                    nodegree + re74 + re75, 
                    data = lalonde,
                    method = "nearest",
                    distance= 'mahalanobis',
                    replace=T,
                    ratio=5
                  )

Once we have our matchit object, the Cobalt package gives us some options for balance checking. After examining the results, you might want to go back and re-specify the model.

bal.tab(ma_matched, thresholds = c(m = .1), un = TRUE)
Balance Measures
               Type Diff.Un Diff.Adj        M.Threshold
age         Contin. -0.3094   0.2904 Not Balanced, >0.1
educ        Contin.  0.0550  -0.0403     Balanced, <0.1
married      Binary -0.3236   0.0314     Balanced, <0.1
race_black   Binary  0.6404   0.0086     Balanced, <0.1
race_hispan  Binary -0.0827   0.0000     Balanced, <0.1
race_white   Binary -0.5577  -0.0086     Balanced, <0.1
nodegree     Binary  0.1114   0.0130     Balanced, <0.1
re74        Contin. -0.7211   0.0817     Balanced, <0.1
re75        Contin. -0.2903   0.1488 Not Balanced, >0.1

Balance tally for mean differences
                   count
Balanced, <0.1         7
Not Balanced, >0.1     2

Variable with the greatest mean difference
 Variable Diff.Adj        M.Threshold
      age   0.2904 Not Balanced, >0.1

Sample sizes
                     Control Treated
All                   429.       185
Matched (ESS)          74.59     185
Matched (Unweighted)  196.       185
Unmatched             233.         0
bal.plot(ma_matched)

love.plot(ma_matched)

Access the matched data with the match_data function:

ma_data<-match_data(ma_matched)

head(ma_data)

Note that a fairly large number of observations have been discarded here:

nrow(ma_data)
[1] 381

We can use the results in our regression model. Be sure to include the weights!

ma_fit<-lm(re78 ~ treat + age + educ +  married + race +nodegree + re74 + re75, 
   data = ma_data,
   weight = weights )

Finally, we want to estimate our effects with robust standard errors:

avg_comparisons(ma_fit,
                variables = "treat",
                vcov = 'HC3'
                )

 Estimate Std. Error    z Pr(>|z|)   S 2.5 % 97.5 %
     1152        900 1.28    0.201 2.3  -612   2915

Term: treat
Type: response
Comparison: 1 - 0
modelsummary(list(
                  "Mahalanobis" = ma_fit
                  ),
             estimate  = "{estimate}",  
             fmt = fmt_significant(),
             statistic ='conf.int',
             conf_level = .95,
             vcov ='HC3',
             gof_omit = 'F|RMSE|R2$|AIC|Log.Lik.', # remove some model stats 
             # add a title
             title = "DV: Wages"
             )
tinytable_urc1y2ir16h3lakobry6
DV: Wages
Mahalanobis
(Intercept) -2890
[-10778, 4998]
treat 1152
[-618, 2921]
age -20.8
[-117.8, 76.2]
educ 703
[117, 1290]
married 160
[-1937, 2257]
racehispan 1718
[-1243, 4678]
racewhite 1246
[-588, 3080]
nodegree 885
[-2036, 3807]
re74 0.0772
[-0.2979, 0.4524]
re75 0.200
[-0.192, 0.593]
Num.Obs. 381
R2 Adj. 0.035
BIC 7982.9
Std.Errors HC3

Coarsened Exact matching

For CEM, we’ll often want to set manually set some of “cutpoints”. There’s no set rules here. Visual inspection of the predictors can help, and the MatchIt package has some defaults that can be sensible, but this is partly a judgement call: what values of the independent variables seem like they should be grouped together?

By default, observations are matched exactly on factor variables. However, you can override this by setting passing a list of levels to the grouping argument:

cem_match <- matchit(treat ~ age + educ +  married + race +
                    nodegree + re74 + re75, data = lalonde,
                  method = "cem",
                  estimand ='ATE',
                  cutpoints =list(educ = c(0,9, 12, 14)),
                  grouping = list(race = list(c("white", "hispan"),
                                              c("black")))
                  )
summary(cem_match)

Call:
matchit(formula = treat ~ age + educ + married + race + nodegree + 
    re74 + re75, data = lalonde, method = "cem", estimand = "ATE", 
    cutpoints = list(educ = c(0, 9, 12, 14)), grouping = list(race = list(c("white", 
        "hispan"), c("black"))))

Summary of Balance for All Data:
           Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
age              25.8162       28.0303         -0.2419     0.4400    0.0813
educ             10.3459       10.2354          0.0448     0.4959    0.0347
married           0.1892        0.5128         -0.7208          .    0.3236
raceblack         0.8432        0.2028          1.6708          .    0.6404
racehispan        0.0595        0.1422         -0.2774          .    0.0827
racewhite         0.0973        0.6550         -1.4080          .    0.5577
nodegree          0.7081        0.5967          0.2355          .    0.1114
re74           2095.5737     5619.2365         -0.5958     0.5181    0.2248
re75           1532.0553     2466.4844         -0.2870     0.9563    0.1342
           eCDF Max
age          0.1577
educ         0.1114
married      0.3236
raceblack    0.6404
racehispan   0.0827
racewhite    0.5577
nodegree     0.1114
re74         0.4470
re75         0.2876

Summary of Balance for Matched Data:
           Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
age              21.1665       20.7539          0.0451     0.8970    0.0144
educ             10.1082       10.1923         -0.0340     0.7899    0.0153
married           0.0619        0.0619         -0.0000          .    0.0000
raceblack         0.5876        0.5876          0.0000          .    0.0000
racehispan        0.1314        0.0940          0.1256          .    0.0375
racewhite         0.2809        0.3184         -0.0946          .    0.0375
nodegree          0.7423        0.7423          0.0000          .    0.0000
re74            472.0782      742.1934         -0.0457     1.2640    0.0478
re75            466.2790      677.2163         -0.0648     1.1178    0.0549
           eCDF Max Std. Pair Dist.
age          0.1294          0.1133
educ         0.1186          0.3455
married      0.0000          0.0000
raceblack    0.0000          0.0000
racehispan   0.0375          0.2567
racewhite    0.0375          0.1933
nodegree     0.0000          0.0000
re74         0.2828          0.0847
re75         0.2699          0.1654

Sample Sizes:
              Control Treated
All            429.    185.  
Matched (ESS)   86.99   32.26
Matched        112.     82.  
Unmatched      317.    103.  
Discarded        0.      0.

Now we can plot our results to check out how well the matching performed:

plot(cem_match, 
     type = "density", 
     interactive = TRUE, 
     which.xs = ~age +
married + re75)

bal.tab(cem_match, thresholds = c(m = .1), un = TRUE)
Balance Measures
               Type Diff.Un Diff.Adj    M.Threshold
age         Contin. -0.2419   0.0451 Balanced, <0.1
educ        Contin.  0.0448  -0.0340 Balanced, <0.1
married      Binary -0.3236  -0.0000 Balanced, <0.1
race_black   Binary  0.6404   0.0000 Balanced, <0.1
race_hispan  Binary -0.0827   0.0375 Balanced, <0.1
race_white   Binary -0.5577  -0.0375 Balanced, <0.1
nodegree     Binary  0.1114   0.0000 Balanced, <0.1
re74        Contin. -0.5958  -0.0457 Balanced, <0.1
re75        Contin. -0.2870  -0.0648 Balanced, <0.1

Balance tally for mean differences
                   count
Balanced, <0.1         9
Not Balanced, >0.1     0

Variable with the greatest mean difference
 Variable Diff.Adj    M.Threshold
     re75  -0.0648 Balanced, <0.1

Sample sizes
                     Control Treated
All                   429.    185.  
Matched (ESS)          86.99   32.26
Matched (Unweighted)  112.     82.  
Unmatched             317.    103.  
bal.plot(cem_match)

love.plot(cem_match)

Finally, we’ll want to use the match_data function to extract the matched data and then estimate our model with the weights included:

cem_data<-match_data(cem_match)


cem_fit <- lm(re78 ~ 
            treat + age + educ + married + re74 +re75+
            nodegree+ race , data = cem_data, weights = weights)
avg_comparisons(cem_fit, variables = "treat", vcov ='HC3', cluster=~subclass)

 Estimate Std. Error    z Pr(>|z|)   S 2.5 % 97.5 %
     1672       1158 1.44    0.149 2.7  -597   3941

Term: treat
Type: response
Comparison: 1 - 0
modelsummary(list(
                   "Mahalanobis" = ma_fit,
                  "Coarsened" = cem_fit
                  ),
             estimate  = "{estimate}",  
             fmt = fmt_significant(),
             statistic ='conf.int',
             conf_level = .95,
             vcov ='HC3',
             gof_omit = 'F|RMSE|R2$|AIC|Log.Lik.', # remove some model stats 
             # add a title
             title = "DV: Wages"
             )
tinytable_anlq9n57jwpcf99n9ghz
DV: Wages
Mahalanobis Coarsened
(Intercept) -2890 -4693
[-10778, 4998] [-16885, 7499]
treat 1152 1672
[-618, 2921] [-612, 3956]
age -20.8 106
[-117.8, 76.2] [-104, 317]
educ 703 608
[117, 1290] [-164, 1379]
married 160 -2533
[-1937, 2257] [-6329, 1263]
racehispan 1718 957
[-1243, 4678] [-1046, 2959]
racewhite 1246 2391
[-588, 3080] [-853, 5636]
nodegree 885 826
[-2036, 3807] [-2735, 4387]
re74 0.0772 -0.0964
[-0.2979, 0.4524] [-1.1358, 0.9430]
re75 0.200 0.687
[-0.192, 0.593] [-0.214, 1.588]
Num.Obs. 381 194
R2 Adj. 0.035 0.061
BIC 7982.9 3999.8
Std.Errors HC3 HC3

Inverse probability weighting

Finally, we can try using inverse probability weighting on the propensity scores instead of matching. Here’s how we might do that using mostly base R functions:

# Step 1: Estimate propensity scores
fit1 <- glm(treat ~age + educ + nodegree + 
                married + race + re74 + re75 , 
            family = binomial, data =  lalonde)
lalonde$propensity <- predict(fit1, type = "response")

# calculate the inverse
lalonde<-lalonde|>
  mutate(ipw = (treat/ propensity) + ((1 - treat) / (1 - propensity)))



# Step 2: Fit weighted outcome model

m <- lm(re78 ~ treat , data = lalonde, weight = ipw)

summary(m)

Call:
lm(formula = re78 ~ treat, data = lalonde, weights = ipw)

Weighted Residuals:
   Min     1Q Median     3Q    Max 
-42083  -6606  -2284   4979  77674 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   6422.8      397.4  16.161   <2e-16 ***
treat          224.7      577.7   0.389    0.697    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9864 on 612 degrees of freedom
Multiple R-squared:  0.0002471, Adjusted R-squared:  -0.001386 
F-statistic: 0.1513 on 1 and 612 DF,  p-value: 0.6975
avg_comparisons(m, variables = "treat", wts = lalonde$ipw, vcov = 'HC0')

 Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
      225        909 0.247    0.805 0.3 -1558   2007

Term: treat
Type: response
Comparison: 1 - 0

Or we can do this using the WeightIt package:

W <- weightit(treat ~ age + educ + nodegree + 
                married + race + re74 + re75, 
              data = lalonde, method = "glm", 
              estimand = "ATE")
fit <- lm_weightit(re78 ~ treat, data = lalonde,
                   weightit = W)

summary(fit, ci = TRUE)

Call:
lm_weightit(formula = re78 ~ treat, data = lalonde, weightit = W)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   2.5 %  97.5 %    
(Intercept)   6422.8      353.4  18.177   <1e-06  5730.3  7115.4 ***
treat          224.7      876.2   0.256    0.798 -1492.6  1942.0    
Standard error: HC0 robust (adjusted for estimation of weights)

Or trying using a different weighting method:

W <- weightit(treat ~ age + educ + nodegree + 
                married + race + re74 + re75, 
              data = lalonde, method = "ebal", 
              estimand = "ATE")
fit <- lm_weightit(re78 ~ treat, data = lalonde,
                   weightit = W)

summary(fit, ci = TRUE)

Call:
lm_weightit(formula = re78 ~ treat, data = lalonde, weightit = W)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)   2.5 %  97.5 %    
(Intercept)   6327.5      326.2  19.400   <1e-06  5688.3  6966.8 ***
treat          951.7     1232.6   0.772     0.44 -1464.2  3367.5    
Standard error: HC0 robust (adjusted for estimation of weights)
avg_comparisons(fit, variables = "treat")

 Estimate Std. Error     z Pr(>|z|)   S 2.5 % 97.5 %
      952       1233 0.772     0.44 1.2 -1464   3368

Term: treat
Type: probs
Comparison: 1 - 0